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FLUKA is a general purpose Monte Carlo transport and interaction code used for fundamental physics and for 
a wide range of applications. These include Cosmic Ray Physics (muons, neutrinos, EAS, underground physics), 
both for basic research and applied studies in space and atmospheric flight dosimetry and radiation damage. A 
review of the hadronic models available in FLUKA and relevant for the description of cosmic ray air showers is 
presented in this paper. Recent updates concerning these models are discussed. The FLUKA capabilities in the 
simulation of the formation and propagation of EM and hadronic showers in the Earth's atmosphere are shown. 



1. Introduction 

Extended Air Showers (EAS) are originated by 
highly energetic cosmic rays, which interact with 
air atoms in the Earth's atmosphere, produc- 
ing elementary hadrons and nuclear fragments. 
Neutral pions immediately decay in two photons, 
which in turn interact with other particles, ge- 
nerating electromagnetic cascades, while charged 
pions and kaons, as well as other hadrons, inter- 
act and/or decay, depending on their energy and 
air density, leading to the hadronic component 
of the shower. Thus, to understand a process 
as complex as an air shower, models and codes 
capable of describing the evolution of both the 
electromagnetic and the hadronic component are 
needed. Monte Carlo codes have been developed 
in the last few years for this purpose, such as 
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the CORSIKA simulation package [T] , allowing to 
choose among different models and model combi- 
nations for the description of EAS formation and 
propagation. 

In this paper, the FLUKA Monte Carlo code [H 
[3] is applied to the study of EAS induced by cos- 
mic rays with primary energies up to 1000 TeV. 
These results rely entirely on the use of the 
FLUKA code for transport, interactions and 
decays of leptons and 7s at all energies and 
for hadron-nucleus interactions at energies below 
20 TeV, while the DPMJET code @E] is used for 
all nucleus-nucleus interactions and for hadron- 
nucleus interactions at the highest energies. Our 
simulations are completely independent with re- 
spect to CORSIKA, which contains only a partial 
version of FLUKA, for use at low energies only (E 
< 100-200 GeV). The energy range spanned in 
this work is of interest for cosmic ray experiments 
aiming at the determination of the primary spec- 
tra and composition for energies up to the knee 



2 



region, such as the ATIC, the RUNJOB and the 
KASCADE experiments. 

It is worth mentioning that FLUKA has not 
been specifically designed for the study of Cosmic 
Ray Physics, but can be applied also in this field. 
At energies < 30 TeV, interesting results con- 
cerning the prediction of inclusive lepton fluxes 
have already been obtained. In particular, the 
first 3-dimcnsional calculation of the atmospheric 
neutrino flux due to Cosmic Rays [5] was made 
with FLUKA, and, more recently, data concer- 
ning muons detected by the balloon-borne BESS 
spectrometer at various depths and by the L3+C 
spectrometer located at CERN have been quite 
succesfully reproduced [7j. 

2. Validation of the FLUKA models 

In this section some aspects of the FLUKA 
models relevant for Cosmic Ray description are 
briefly discussed, together with the validation of 
the models by means of data collected at acce- 
lerators. Actually, even the most recent accele- 
rator data allow to test theoretical Monte Carlo 
models only for some specific projectile beam - 
target combinations at energies well below the 
maximum primary energies, while, for the simu- 
lation of EAS processes, one has to extrapolate 
the models to many more combinations and to 
the highest energies, on the basis of theoretical 
considerations. 

2.1. FLUKA: a brief introduction 

The modern FLUKA |2|3j is a general purpose 
Monte Carlo transport and interaction code de- 
veloped since 1988 mostly by INFN and CERN 
researchers, and used in many fields of physics, 
both for fundamental research and for applica- 
tions. The history of the code dates back to the 
sixties when J. Ranft developed and applied the 
very first version to accelerator shielding. The 
most common FLUKA applications are accelera- 
tor physics, calorimetry, shielding design, radio- 
protection and dosimetry, while new ones, such as 
hadrontherapy, are promisingly growing. Since a 
few years, FLUKA is also applied to Cosmic Ray 
Physics, mainly with focus on low energy cosmic 
rays, as already mentioned in the Introduction, 
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Figure 1. Computed ir + double differential pro- 
duction cross section for 12.9 GeV/c protons on 
Aluminum for different angular ranges, compared 
with HARP experimental data [19]. 



and on their effects for dosimetry |8I9| and space 
physics [TO]- All these applications require both 
precise physics models, and packages to build ge- 
ometries, to simulate targets and experimental 
setups. Both topics are the subject of constant 
development and improvement in FLUKA. 

As far as the geometry is concerned, a combina- 
torial geometry package is included into the code, 
allowing to create complex geometries, which 
can be visualized by means of apposite graphical 
tools. This allows to simulate 3-D particle propa- 
gation, and, as far as Cosmic Rays are concerned, 
to study showers characterized by whichever in- 
clination angle and direction with respect to the 
Earth, and to take into account in a precise way 
the effect of different atmospheric compositions 
and magnetic field configurations. The same code 
can be used for the simulation of cosmic ray de- 
tectors, including surrounding materials. 

The FLUKA physical models are based as far 
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Figure 2. Feynman X distributions for 7r + and 
7r~ production for proton interactions on Carbon 
at 158 GeV/c , as measured by NA49 [20] (sym- 
bols) and predicted by FLUKA (histograms). 
Linear scale on the left, logarithmic scale on the 
right. 



cial treatment of the tip of the bremsstrahlung 
spectrum. Electron pairs and bremsstrahlung 
are sampled from the proper double differential 
energy-angular distributions improving the com- 
mon practice of using average angles. In a similar 
way, the three-dimensional shape of the hadronic 
cascades is reproduced in detail by a rigorous 
sampling of correlated energy and angles in de- 
cay, scattering, and multiple Coulomb scattering. 

Bremsstrahlung and direct pair production by 
muons arc modelled according to state-of-the- 
art theoretical description and have been checked 
against experimental data |12|13j . Muon pho- 
tonuclear interactions are also modelled. 

2.3. The FLUKA hadronic models 



as possible on theoretical microscopic models, 
with the advantage, with respect to paramete- 
rized inclusive models, of preserving correlations, 
and of being predictive in the regions where ex- 
perimental data are not available. The model pa- 
rameters are fixed once, internally to the code, for 
all projectile-target combinations and energies, 
and cannot be modified by the user. The theoreti- 
cal models are continuously benchmarkcd against 
newly available experimental data. The code is 
under continuous development and made availa- 
ble on the website http://www.fluka.org . The de- 
velopment and maintenance are performed under 
a INFN-CERN agreement. 

2.2. E.M. and muon transport in FLUKA 

For historical reasons, FLUKA is best known 
for its hadron event generators, but since more 
than 17 years FLUKA can handle with similar 
or better accuracy electromagnetic effects [TT] . 
Briefly, the energy range covered by this sector 
of FLUKA is very wide: the program can trans- 
port photons and electrons over about 12 energy 
decades, from 1 PeV down to 1 keV. The e.m. 
part is fully coupled with the hadron sector, in- 
cluding the low energy (i.e. < 20 MeV) neu- 
trons. The simulation of the electromagnetic ca- 
scade in FLUKA is very accurate, including the 
Landau-Pomeranchuk-Migdal effect and a spe- 




Figure 3. Average e/fi fluence as a function of the 
atmospheric depth for vertical showers induced by 
p, Fe primaries with 10 14 eV (upper panel) and 
10 15 eV (lower panel) energies. In each figure e 
from Fe (asterisks), e fromp (+ symbols), fi from 
Fe (filled squares) and \i from p (filled circles) are 
shown. The plots refer to e in the energy range 1 
MeV - 1 TeV and to /i in the energy range 1 GeV 
- 1 TeV. 
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Figure 4. /i fluence as a function of the at- 
mospheric depth for vertical showers induced by 
p (upper panel) and Fe (lower panel) primaries 
with 10 15 eV energy. Each line in each figure 
is the result for a different shower. In case of p 
showers, the development of the hadronic compo- 
nent shows larger fluctuations, also related to the 
depth of the first interaction. See also Fig. [31 



A basic description of hadronic interactions in 
FLUKA can be found in [14] . Hadron-nucleon 
interactions at energies below a few GeV are si- 
mulated in FLUKA by the isobar model, through 
resonance production and decay, and by taking 
into account elastic, charge and strangeness ex- 
change. Elementary hadron-hadron collisions at 
energies above a few GeV are described thanks 
to an implementation of the Dual Parton Model 
(DPM) [15], coupled to a hadronization scheme. 
The Dual Parton Model allows for the description 
of soft collision processes (i.e. processes characte- 
rized by low pt,Pt « Qo, where Qo is a momen- 
tum scale of the order of a few GeV) , that cannot 
be described by means of pQCD, due to fast rise 
of as with decreasing momentum. DPM incor- 
porates analiticity, duality and unitarity. Incom- 
ing hadrons are described by strings, which inter- 
act by exchanging closed loop excitations, called 



Figure 5. Spatial distribution of e fluence 
(particles/cm 2 /primary) for vertical showers in- 
duced by Fe primaries with 10 13 eV (upper panel) 
and 10 15 eV (lower panel) energies. The results 
are obtained as an average on ~ 100 events for 
each energy. The primaries come from the top of 
the atmosphere (top of each figure) and propaga- 
te towards the Earth's surface located at ~ 6378 
km with respect to the Earth's center. 



pomerons. At low energies, the single-pomeron 
exchange is the dominant contribution, while at 
laboratory energies >> 1 TeV higher-order con- 
tributions, corresponding to multi-pomeron ex- 
change, become increasingly important. Corre- 
sponding to each pomeron cut, built according to 
the recipes provided by the Abramoviskii-Gribov- 
Kancheli cutting rules, colorless chains are built, 
extending from the projectile valence and sea 
quarks and antiquarks to the target ones. In 
particular, the valence quarks in each baryon 
are treated as a quark-diquark pair, so that, 
in the simplest case of the single-pomeron ex- 
change amplitude for a baryon-baryon interac- 
tion, two chains are built, each extending from 
a valence quark to a valence diquark. Sea quarks 
(and the corresponding antiquarks) are additio- 
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Figure 6. Same as Fig. [5] for fi fluence 
(particles/cm 2 /primary) in vertical showers in- 
duced by Fe primaries with 10 13 eV (upper panel) 
and 10 15 eV (lower panel) energies. 



Figure 7. Spatial distribution of EM energy de- 
position (GeV/cm 3 /primary) for vertical showers 
induced by Fe ion primaries with 10 13 eV (upper 
panel) and 10 15 eV (lower panel) incoming ener- 
gies. See also Fig. [5] 



nally involved as end-points of chains in multi- 
pomeron exchange amplitudes. Besides determi- 
ning the number of pomeron cuts and forming 
the chains, the DPM allows to assign an energy 
and momentum to each chain, according to the 
momentum distribution functions of the partons 
in the hadrons. In the asymptotic regime par- 
ton masses are neglected, however, at lower ener- 
gies, finite mass effects have to be included both 
in the chain building process, by suppressing the 
chains with invariant masses below the observed 
baryon and meson mass values and reassigning 
energy and momentum to other chains to en- 
sure energy/momentum conservation, and in the 
hadronization procedure. Hadronization is not 
included in DPM, but is performed by different 
models properly coupled to DPM. At present, the 
hadronization scheme implemented in FLUKA is 
based on an advanced version [TH] of the model 



by [H] , with special emphasis to provide accuracy 
in the low-energy description of hadron formation 
(down to a few GeV in the lab). 

Hadron-hadron collisions are the main build- 
ing blocks of hadron- nucleus collisions. Multiple 
collisions of each hadron with the nuclear con- 
stituents are taken into account by means of the 
Glauber-Gribov mechanism. Particular efforts 
are devoted to the study of nuclear effects on 
hadron propagation. In the last FLUKA version 
the propagation in the nucleus of the hadrons re- 
sulting from elementary multiple collisions can be 
simulated with improved accuracy with respect 
to the past, by means of an improved Genera- 
lized IntraNuclear Cascade (GINC) model, fol- 
lowed by pre-equilibrium and coalescence and by 
de-excitation of the remaining nuclear fragments. 
All these models are included in the FLUKA sub- 
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module PEANUT, which has been refined for 
several years to describe with increasing accu- 
racy many low-energy processes involving nucle- 
ons, and has been further extended in the last 
few months to higher energy events (E > a few 
GeV). More detail about PEANUT and the issues 
concerning its extension can be found in [14|18j . 
PEANUT is used also to simulate 7, v and stop- 
ping [i interactions, and nucleon decay. Photon 
interactions with hadrons and muon virtual pho- 
ton nuclear interactions are simulated by means 
of the Vector Dominance Model. 

Hadron free paths inside a nucleus can be dif- 
ferent from those computed on the basis of free 
hadron-nucleon scattering, due to Pauli blocking, 
coherence length, formation and multibody pro- 
cesses, all implemented in FLUKA |14|18j . In 
particular, coherence length affects final states 
from elastic and charge-exchange scatterings: in- 
teraction products can not be localized better 
than the position uncertainty related to the 4- 
momentum transferred in the collision, according 
to the AxAp > h relationship. Thus, reinterac- 
tions occurring at distances shorter than the co- 
herence length undergo interference and can not 
be treated as independent. On the other hand, 
formation zone affects the reinteraction probabi- 
lity of secondaries emerging from high energy in- 
teractions. It is supported by experimental re- 
sults, which indicate that high-energy secondary 
particle reinteractions inside nuclei are strongly 
suppressed. Taking into account that a typical 
time for strong interactions is ~ 1 fm/c, a par- 
ticle emerging from an interaction requires some 
time to materialize, according to the uncertainty 
principle. 

Examples of performance of the FLUKA 
hadronic model, in reproducing experimental 
data recently obtained by the HARP and NA49 
Collaborations are shown in Fig. Q] and Fig. [2] 

2.4. The FLUKA-DPMJET interface 

Both hadronic physics at energies above 20 TeV 
and heavy-ion collisions above 5 GeV have been 
simulated in this work thanks to the interface 
of FLUKA with the DPMJET code. DPMJET- 
11.53 and DPMJET-III libraries can be both used 
within FLUKA. DPMJET is also available in 
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Figure 8. (N e , A M ) correlation at the Earth sea 
level: dependence on the primary mass and en- 
ergy. Each symbol corresponds to a different 
shower event. Blobs from the lower left corner to 
the upper right corner refer to events originated 
from vertical p and Fe ions with initial energies 
10 12 , 10 13 , 10 14 and 10 15 eV, respectively. 
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Figure 9. Same as Fig. [8] at a 4000 m a.s.l. alti- 
tude level. Calculations at such high altitudes are 
of interest for EAS experiments located at moun- 
tain sites. 



CORSIKA as a possible choice among multiple 
high-energy event generators, and is widely used 
in the Cosmic Ray Physics community. DPMJET 
has been extensively tested against data from the 
SPS collider, the Tevatron and RHIC. 

Soft physics is described in DPMJET gE] 
thanks to the DPM, already introduced in the 
previous section. Furthermore, hard physics pro- 
cesses are implemented, based on leading order 
pQCD. Soft and hard processes contribute to- 
gether to the inelastic cross-section according to 
the cikonal approximation, which provides a uni- 
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tarization scheme. The eikonal function is given 
by the sum of the soft and hard components, also 
including terms from multi-pomeron exchange. 

One of the main differences between the 
FLUKA hadronic physics models and DPMJET 
is related to the hadronization process. In DPM- 
JET this process is performed by means of the 
JETSET chain fragmentation code included in 
PYTHIA [21] , with some modifications. In par- 
ticular, while the PYTHIA code was extensively 
tested for processes involving a significant hard 
component, the transverse momentum distribu- 
tion in soft chain decay has been modified in the 
JETSET version included in DPMJET to bet- 
ter describe soft chain hadronization. This way, 
the transverse momentum distributions of K and 
p predicted by DPMJET-III show an improved 
agreement with experimental data. 

There are other differences in the description 
of final state interactions: DPMJET has its 
own GINC code, simplified with respect to the 
FLUKA one, and does not take into account 
coalescence and pre-equilibrium emission. Nu- 
clear de-excitation is performed here by means 
of the same de-excitation models [14122123] used 
in FLUKA. 

Further refinements of DPMJET and of its in- 
terface to FLUKA are under way. As an example, 
diquark breaking mechanisms, absent in the ori- 
ginal DPM, have recently been proposed, in the 
attempt of providing a contribution to the en- 
hanced stopping experimentally observed for nu- 
clear collisions in fixed target experiments. The 
effects on the prediction of particle/antiparticle 
asymmetries, which can be compared with mea- 
sured ones, are under investigation. 

3. Theoretical EAS predictions 

To test the capabilities of our model, we have 
simulated hundreds of vertical air shower events, 
for primary protons and iron nuclei, with ener- 
gies in the 10 12 — 10 15 eV range. The average 
longitudinal EAS e//x profiles are shown in the 
upper and lower panels of Fig.[3]for 10 14 and 10 15 
eV primary energies, respectively. Electrons and 
positrons of energies > 1 MeV have been consid- 
ered, together with muons of energies > 1 GeV, 



up to 1 TeV. As can be seen from each figure, 
at a fixed energy, the X^ux depth decreases with 
increasing primary mass, due to the behaviour of 
the p-air inelastic scattering cross-section, which 
indeed increases with mass number. The same 
is true also for /i^ , whose longitudinal profile is 
however softer, i.e. has a maximum less peaked 
than the one. More detail about the /i longi- 
tudinal development can be inferred from Fig. 2] 
where profiles relative to individual showers in- 
duced by p (upper panel) and Fe (lower panel) 
primaries at 10 15 eV are shown. Profile fluctua- 
tions are more evident for lower mass primaries, 
also related to the location of the first primary in- 
teraction with air, which can even vary within a 
range of a few km for low mass primaries and low 
density atmosphere. The muon number reflects 
the development of the hadronic component of 
the shower, and is thus a basic test of the hadronic 
models, since muons come from charged pion and 
kaon decay. These processes become important at 
energies below the critical energy, i.e. the energy 
for which the meson interaction and decay lengths 
are nearly the same. In particular the critical en- 
ergy for 7r decay is around 150 GeV, while the 
one for K is about ten times larger. From these 
considerations, it is clear that the /i^ number is 
particularly sensitive to a good description of the 
hadronic physics processes at low energies. We 
emphasize that we refer indeed to with en- 
ergy down to 1 GeV. On the other hand, with 
energies around hundreds or thousands GeV in 
the atmosphere deeply penetrate the rock, and 
can be selected and detected by underground de- 
tectors. Multi-// and /^-bundle events have also 
been studied with FLUKA, which has provided a 
succesful reproduction of experimental data [7] . 

To better identify the longitudinal shape of the 
shower components, detailed maps of fluences and 
energy deposition have been built. In Fig. [5] and [6] 
the fluences of and /i^ from the top of the at- 
mosphere, down to the Earth's surface, are shown 
for Fe induced vertical showers at 10 13 and 10 15 
eV. From the figures, it is evident that in the top 
part of the atmosphere the e^- shower profile is 
more extended than the ^ one, and that some 
can even propagate from lower altitudes to- 
wards higher ones. At the considered energies, 
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some and ^ from the same shower can reach 
ground within a few km. Anyway, the largest flu- 
ences are found for both kind of particles in the 
spatial region located within ^2-8 km above 
the Earth's surface, at distances not larger than 
500 - 700 m from the shower axis. In Fig. [7] de- 
tailed maps of electromagnetic energy deposition 
for vertical Fe showers at 10 13 and 10 15 eV are 
presented. These maps are strictly related to the 
ones for e fluence shown in Fig. [5l already discus- 
sed. e ± production thresholds have been set to 1 
MeV, while 7 threshold has been set to 0.5 MeV 
in our simulations. 

The (N e , N^) correlations for p and Fe show- 
ers at sea level and at 4000 m a.s.l. are shown 
in Fig. [8] and Fig. [9l respectively. Results for pri- 
mary energies from 10 12 to 10 15 eV are plotted. It 
appears that for lower energy showers and higher 
altitudes it is more difficult to distinguish the pri- 
mary spectrum composition on the basis of (N e , 
N^) detection. 

The features of EAS predicted by our simula- 
tions are in a qualitative agreement with those 
obtained by other authors. As a further step in 
our investigation, we plan to compare in more 
detail our theoretical results with those obtained 
by means of other models, e.g. those available in 
the CORSIKA code. Beyond the question of the 
comparison of hadronic interaction models, this 
can be helpful to identify also possible differences 
coming from different e.m. or transport models. 
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